
rm(list = ls()) #CLEARS GLOBAL ENVIRONMENT
cat ("\014")    #CLEARS THE R CONSULE


#######################################
## Agency cooperation and reputation ##
#######################################

# to read CSV files
#install.packages("lavaan", dependencies = TRUE)
library(lavaan)
#install.packages("foreign")
library(foreign)
#install.packages("dplyr")
library(dplyr)

#install.packages("lmtest")
library(lmtest)
#install.packages("sandwich")
library(sandwich)
#install.packages('modelsummary')
library(modelsummary)
library(car)
#install.packages("plm")
library(plm)
#install.packages("estimatr")
library(estimatr)
#install.packages("texreg")
library(texreg)
#install.packages("xtable")
library(xtable)
library(tidyverse)
#install.packages("stargazer")
library(stargazer)

#READ DATA
coop_rep <- read.csv("C:/[set path]/Colab_Rep.csv")

#DELETE IRRLEVANT AGENCIES
coop_rep <- filter(coop_rep, name!="SIPC")
coop_rep <- filter(coop_rep, name!="NIBS")

#CREATE REPUTATION VALUE (NATURAL LOG OF POS-NEG RATIO):
coop_rep <- transform(coop_rep,
            rep = log((positive/negative)+1)
)

#=========================================================
# REPUTATION BY COOPERATION

# AUTOCORRELATION MODELS:
AR_rep <- plm(rep ~
                lag(rep),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(AR_rep)

AR_colab <- plm(d_comulative_MOU ~
                lag(d_comulative_MOU),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(AR_colab)

#WITHIN ESTIMATOR (FE MODEL):
model1 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU)+
                lag(rep,2),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model1)

model2 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,2)+
                lag(rep,3),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model2)

model3 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,3)+
                lag(rep,4),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model3)

model4 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,4)+
                lag(rep,5),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model4)

rob_se <- list(sqrt(diag(vcovHC(model1, type = "HC1"))),
               sqrt(diag(vcovHC(model2, type = "HC1"))),
               sqrt(diag(vcovHC(model3, type = "HC1"))),
               sqrt(diag(vcovHC(model4, type = "HC1"))))

stargazer(model1, model2, model3, model4, title='Results',
          type="text",
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          style = "ajps",
          se = rob_se,
          out = "table_XA.html",
          keep=c("d_comulative_MOU","rep"))

#cOEFFICIENT PLOT:
models <- list()
models[['1y lagged collaboration']] <- model1
models[['2y lagged collaboration']] <- model2
models[['3y lagged collaboration']] <- model3
models[['4y lagged collaboration']] <- model4

modelplot(models,
          coef_omit = 'rep') + theme_classic()

#===============================================
# COOPERATION BY REPUPTION:

model5 <- plm(d_comulative_MOU ~
                lag(rep)+
                lag(d_comulative_MOU,2)+
                rep,
              index="numagency",
              model="within",
              data = coop_rep)
summary(model5)

model6 <- plm(d_comulative_MOU ~
                lag(rep,2)+
                lag(d_comulative_MOU,3)+
                rep,
              index="numagency",
              model="within",
              data = coop_rep)
summary(model6)

model7 <- plm(d_comulative_MOU ~
                lag(rep,3)+
                lag(d_comulative_MOU,4)+
                rep,
              index="numagency",
              model="within",
              data = coop_rep)
summary(model7)

model8 <- plm(d_comulative_MOU ~
                lag(rep,4)+
                lag(d_comulative_MOU,5)+
                rep,
              index="numagency",
              model="within",
              data = coop_rep)
summary(model8)

rob_se_r <- list(sqrt(diag(vcovHC(model5, type = "HC1"))),
                 sqrt(diag(vcovHC(model6, type = "HC1"))),
                 sqrt(diag(vcovHC(model7, type = "HC1"))),
                 sqrt(diag(vcovHC(model8, type = "HC1"))))

stargazer(model5, model6, model7, model8, title='Results',
          type="text",
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          style = "ajps",
          se = rob_se_r,
          out = "table_XB.html",
          keep=c("d_comulative_MOU","rep"))


#===============================================
#MODELS 1-4 WITH AGENCIES THAT HAD AT LEAST SOME COLLABORATION WITH THE FDA & controlling for % articles that include the FDA - TABLE A1

#DELETING AGENCIES WITH NO MOU:
coop_rep <- filter(coop_rep, name!="BPD")
coop_rep <- filter(coop_rep, name!="CFTC")
coop_rep <- filter(coop_rep, name!="DOL")
coop_rep <- filter(coop_rep, name!="ETA")
coop_rep <- filter(coop_rep, name!="FHWA")
coop_rep <- filter(coop_rep, name!="FHA")
coop_rep <- filter(coop_rep, name!="FAS")
coop_rep <- filter(coop_rep, name!="GSA")
coop_rep <- filter(coop_rep, name!="NMFS")
coop_rep <- filter(coop_rep, name!="NTSB")
coop_rep <- filter(coop_rep, name!="OPM")
coop_rep <- filter(coop_rep, name!="OSC")
coop_rep <- filter(coop_rep, name!="SBA")
coop_rep <- filter(coop_rep, name!="USITC")

# 

#WITHIN ESTIMATOR (FE MODEL):
model1 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU)+
                lag(rep,2),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model1)

model2 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,2)+
                lag(rep,3),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model2)

model3 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,3)+
                lag(rep,4),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model3)

model4 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,4)+
                lag(rep,5),
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model4)

rob_se <- list(sqrt(diag(vcovHC(model1, type = "HC1"))),
               sqrt(diag(vcovHC(model2, type = "HC1"))),
               sqrt(diag(vcovHC(model3, type = "HC1"))),
               sqrt(diag(vcovHC(model4, type = "HC1"))))

stargazer(model1, model2, model3, model4, title='Results',
          type="text",
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          style = "ajps",
          se = rob_se,
          out = "table_XA.html",
          keep=c("d_comulative_MOU","rep", "with_FDA"))

# CONTROLLING FOR THE SHARE OF ARTICLES IN WHICH THE FDA APEARS WITH THE COLLABORATING AGENCY - TABLE A2

#WITHIN ESTIMATOR (FE MODEL):
model1 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU)+
                lag(rep,2)+with_FDA,
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model1)

model2 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,2)+
                lag(rep,3)+with_FDA,
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model2)

model3 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,3)+
                lag(rep,4)+with_FDA,
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model3)

model4 <- plm(rep ~
                d_comulative_MOU+
                lag(d_comulative_MOU,4)+
                lag(rep,5)+with_FDA,
              index=c("numagency","year"),
              model="within",
              data = coop_rep)
summary(model4)

rob_se <- list(sqrt(diag(vcovHC(model1, type = "HC1"))),
               sqrt(diag(vcovHC(model2, type = "HC1"))),
               sqrt(diag(vcovHC(model3, type = "HC1"))),
               sqrt(diag(vcovHC(model4, type = "HC1"))))

stargazer(model1, model2, model3, model4, title='Results',
          type="text",
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          style = "ajps",
          se = rob_se,
          out = "table_XA.html",
          keep=c("d_comulative_MOU","rep", "with_FDA"))


#PLOT ESTIMATES
#install.packages("dotwhisker")
library(dotwhisker)
library("ggplot2")
#install.packages("coefplot")
library(coefplot)
#nstall.packages("jtools")
library(jtools)

dwplot(model1, model2, model3, model4)


#REP HISTOGRAM:
library("ggplot2")
ggplot(data = coop_rep, aes(x = rep))+
  geom_histogram()
ggplot(data = coop_rep, aes(x = rep))+
  geom_histogram(binwidth = .1, aes(y = ..density..))+
  geom_density(size = 2, colour = "red")

ggplot(data = coop_rep, aes(x = d_comulative_MOU))+
  geom_histogram()
